#include <iostream>
#include <netcdf.h>
#include "getVar.hpp"
#include <cstring>
#include <math.h>
#include "gcr.h"
#include "norm.hpp"
#include <mkl.h>
#include <mkl_cblas.h>
using namespace std ;
#define ERR {if(status!=NC_NOERR)printf("Error at line=%d: %s %s\n", __LINE__, varname, nc_strerror(status));}
#define ERR1 {if(status!=NC_NOERR)printf("Error at line=%d:  %s\n", __LINE__, nc_strerror(status));}
#define ISNETCDF4ATT "_IsNetcdf4"

void check(int line, int status){
      if (status != NC_NOERR) {
          std::cout << "line "<<line<<", error code " <<status<<": "<<nc_strerror(status) << std::endl;
      }
}
//int data1[5999999], fill_value_in[99999];
void mygetVar_int(int ncid,const char* varname, int * data)
{
    int varid, status;
    check(23, nc_inq_varid(ncid, varname, &varid));
    check(38,nc_get_var_int(ncid,varid,data));
}

extern "C" void module_gcr_mp_solve_helmholts_(float *a_helm,float *b_helm,float *cm_helm,float *threshold, double *pi,  
        int *idep, int *jdep,int *ids, int *ide, int *jds, int *jde, int *kds, int *kde, int *
        ims, int *ime, int *jms, int *jme, int *kms, int *kme, int *
        its, int *ite, int *jts, int *jte, int *kts, int *kte, int *myprcid);
extern "C" void cu_solve_helmholts_(float *a_helm,float *b_helm,float *cm_helm,float *threshold, double *pi, 
        int *idep, int *jdep,int *ids, int *ide, int *jds, int *jde, int *kds, int *kde, int *
        ims, int *ime, int *jms, int *jme, int *kms, int *kme, int *
        its, int *ite, int *jts, int *jte, int *kts, int *kte, int *myprcid, int iter_max, int flag);


extern "C" void module_halo_mp_getijk_(int *);

extern "C" void mygetvar_c_(char postname[9], int *mytid, int *iter_max){
    //string filename="simple_xy.nc";
    int ncid,ncid2, dim_len, ijk_index[20];
    int  idep, jdep, ids,ide,jds,jde,kds,kde,   
                       ims,ime,jms,jme,kms,kme,   
                       its,ite,jts,jte,kts,kte;
    std::string dirname= ("../data/")+string(postname) ;
    //std::cout<<dirname<<dirname.length()<<std::endl;
    check(75,nc_open(dirname.c_str(), NC_NOWRITE, &ncid));
    mygetVar_int(ncid, "ijk_index",ijk_index);

    idep=ijk_index[0]; jdep=ijk_index[1]; ids=ijk_index[2];ide=ijk_index[3];jds=ijk_index[4];jde=ijk_index[5];kds=ijk_index[6];kde=ijk_index[7];
    ims=ijk_index[8];ime=ijk_index[9];jms=ijk_index[10];jme=ijk_index[11];kms=ijk_index[12];kme=ijk_index[13];
    its=ijk_index[14];ite=ijk_index[15];jts=ijk_index[16];jte=ijk_index[17];kts=ijk_index[18];kte=ijk_index[19];
        //std::cout<<"tid:"<<mytid[0]<<" ims-its:"<<ims-its<<", jms-jts:"<<jms-jts<<" ime-ite:"<<ime-ite<<", jme-jte:"<< jme-jte<<std::endl;

    float *a_helm=new float[19*(ite-its+1)*(kte-kts+3)*(jte-jts+1)];
    float *b_helm=new float[(ite-its+1)*(kte-kts+3)*(jte-jts+1)];
    float *b_helm2=new float[(ite-its+1)*(kte-kts+3)*(jte-jts+1)];
    float *cm_helm=new float[7*(ite-its+3)*(kte-kts+3)*(jte-jts+3)];
    float *p1=new float[(ite-its+3)*(kte-kts+3)*(jte-jts+3)];
    float *pi=new float[(ime-ims+1)*(kte-kts+3)*(jme-jms+1)];
//    float *pi2=new float[(ime-ims+1)*(kte-kts+3)*(jme-jms+1)];
    double *pi_r8=new double[(ime-ims+1)*(kte-kts+3)*(jme-jms+1)];
    double *dbl_pi=new double[(ime-ims+1)*(kte-kts+3)*(jme-jms+1)];
//    double *pi_r8_diff=new double[(ime-ims+1)*(kte-kts+3)*(jme-jms+1)];
    mygetVar(ncid, "a_helm",a_helm);
    mygetVar(ncid, "b_helm",b_helm);
    mygetVar(ncid, "cm_helm",cm_helm);
    mygetVar(ncid, "pi",pi);
    mygetVar(ncid, "dbl_pi",dbl_pi);
//    dirname= ("../data/datak/")+string(postname) ;
//    check(75,nc_open(dirname.c_str(), NC_NOWRITE, &ncid2));
//    mygetVar(ncid2, "pi",pi2);
//    mygetVar(ncid2, "b_helm",b_helm2);
    //mygetVar(ncid, "dbl_pi",dbl_pi);
    //for(int i=0;i<7*(ite-its+3)*(kte-kts+3)*(jte-jts+3);i++){
    //  if(fabs(cm_helm[i])<1e-8) printf("%d,%d,",mytid[0],i);
    //}
         module_halo_mp_getijk_(ijk_index);
         //*iter_max=10;
    float threshold=1e-11;
    int i,j,k;
        double c12_tmp;
        int NX = ite - its + 3;
        int NY = kte - kts + 3;
        int NZ = jte - jts + 3;
        int NG1 = NX * NY * NZ;
        int NG2= (ime - ims + 1 ) * NY * (jme - jms + 1);
//#define fortran
/*
#ifdef fortran
    for(i=0;i<(ime-ims+1)*(kte-kts+3)*(jme-jms+1);i++){
        pi_r8[i]=(double)pi[i];
        pi_r8_diff[i]=0;
    }
    //for(int t=0;t<t_max;t++){
    module_gcr_mp_solve_helmholts_(a_helm,b_helm,cm_helm,&threshold, pi_r8,
        &idep, &jdep,&ids, &ide, &jds, &jde, &kds, &kde, 
        &ims, &ime, &jms, &jme, &kms, &kme, 
        &its, &ite, &jts, &jte, &kts, &kte, mytid);
        //if(mytid[0]==12) {
                boundary1(pi_r8_diff[index_x(i,k,j)]=0;)
                para_jte(pi_r8_diff[index_x(i,k,j)]=pi_r8[index_x(i,k,j)]-(dbl_pi[index_x(i,k,j)]);)
                c12_tmp=cblas_ddot(NG2,pi_r8_diff,1,pi_r8_diff,1);
                printf("mytid=%d,diff_norm=%E\n",mytid[0],c12_tmp);
        //}
#endif
*/
//else
//return;
   for(i=0;i<(ime-ims+1)*(kte-kts+3)*(jme-jms+1);i++){ pi_r8[i]=(double)pi[i]; }
   cu_solve_helmholts_(a_helm,b_helm,cm_helm,&threshold, pi_r8, 
        &idep, &jdep,&ids, &ide, &jds, &jde, &kds, &kde, 
        &ims, &ime, &jms, &jme, &kms, &kme, 
        &its, &ite, &jts, &jte, &kts, &kte, mytid, *iter_max, 0);
  return; 
  /* for(i=0;i<(ime-ims+1)*(kte-kts+3)*(jme-jms+1);i++){ pi_r8[i]=(double)pi2[i]; }
   cu_solve_helmholts_(a_helm,b_helm2,cm_helm,&threshold, pi_r8, 
        &idep, &jdep,&ids, &ide, &jds, &jde, &kds, &kde, 
        &ims, &ime, &jms, &jme, &kms, &kme, 
        &its, &ite, &jts, &jte, &kts, &kte, mytid, *iter_max, 1);
    */    
        //if(mytid[0]==12) {
        /*
                boundary1(pi_r8_diff[index_x(i,k,j)]=0;)
                para_jte(pi_r8_diff[index_x(i,k,j)]=pi_r8[index_x(i,k,j)]-(dbl_pi[index_x(i,k,j)]);)
                c12_tmp=cblas_ddot(NG2,pi_r8_diff,1,pi_r8_diff,1);
                if(mytid[0]==5)printf("mytid=%d,diff_norm2=%E\n",mytid[0],c12_tmp);
          */      
        //}
}
